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Abstract Speciation has never been directly observed 
in nature because it is a lengthy phenomenon. 
Although rare, ring species are an optimal natural 
example of speciation and can be identified through 
the assessment of the geographical conditions of their 
potential habitat. Phrynocephalus forsythii is endemic 
to the Tarim Basin, which comprises the Taklamakan 
Desert and surrounded by mountains on three sides. 
This study aimed to determine whether P. forsythii 
had a ring-species-like divergence pattern through the 
characterization of the genetic features of 17 populations 
covering the major distribution of this species. Species 
distribution modelling revealed that P. forsythii had a 
continuous circular distribution around the Tarim Basin. 
Gene flow was observed in most adjacent populations 
except for two terminal populations of the ring, which 
exhibit the highest differentiation. Genetic distance 
and geographic distance were significantly correlated, 
indicating that the observed differentiation resulted 
from genetic variation gradually accumulating during 
population dispersion. Although our results do not 
definitively indicate that P. forsythii is a ring species, 
our results indicate a ring-shaped diversification. This 
phenomenon elucidates the potential mechanism 
underlying speciation in the presence of gene flow, 
providing insight into this evolutionary process. 
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1. Introduction 


Ring species are among the best natural illustrations of 
speciation (Kuchta, 2016). Understanding speciation is a 
fundamental biological problem. During classical ecological 
speciation, population fragmentation yields new species from 
geographically isolated populations of the same ancestral 
species (Dieckmann and Doebeli, 1999; Schluter, 2001). However, 
speciation is gradual phenomenon and thus has not been 
observed directly in nature. Instead, assessment of ring species 
may provide insights into speciation, a chain of intergrading 
populations that encircle a barrier, with terminal forms 
coexisting without interbreeding (Mayr, 1942). Ring species 
represent a phenomenon wherein reproductive isolation can be 
produced in the presence of continuous gene flow and genetic 
variation (Irwin and Wake, 2016; Hey et al, 2005). In brief, ring 
species help us to understand how intraspecific geographical 
variation generates species-level divergence and are integral to 
speciation studies (Irwin et al, 2001; Cacho and Baum, 2012). 
Considering their evolutionary significance, numerous 
studies have attempted to identify ring species. Ideally, the 
criteria that should be met are: circular distribution, gene flow, 
and reproductive isolation (Joseph et al, 2008; Qiao et al., 2018). 
A typical ring species must first meet the requirements in terms 
of distribution, in that the sub-populations around a barrier 
are derived from an ancestral population. Such barriers can 
consist of regions with a slope around them or environmental 
variations that do not generate a suitable habitat for a 
specific species (Irwin, 2012; Fuchs et al., 2015; Pereira and 
Wake, 2015). Nevertheless, gene flow may link the circular 
chain of populations around barriers. With an increase in the 
geographical distance around the ring, genetic variation among 
populations would also gradually accumulate (Alcaide et al, 
2014). This accumulation of genetic differences eventually leads 
to reproductive isolation between the two sub-populations at 
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the far ends of the chain (terminal populations) (Dufresnes et al., 
2016; Irwin et al, 2001). 

A promising candidate for ring species is Phrynocephalus 
forsythii, a lizard endemic to the Tarim Basin (Adler and Zhao, 
1993). Surrounded by mountains on three sides, the Tarim 
Basin also contains the largest shifting-sand desert in central 
Asia, the Taklamakan Desert (Zheng et al., 2009). Furthermore, 
the outer circle of the basin is surrounded by an oasis. These 
environmental conditions provide an ideal environment for the 
evolution of ring-like species. P. forsythii is distributed in the 
oasis around the Tarim Basin (Figure 1) and prefers a habitat 
with dwarf shrubs (Zhang et al., 2010) a humid environment, 
and a lower ambient temperature (Qi et al., 2019). Previous 
studies on the phylogeographic structure of P. forsythii have 
reported that two adjacent populations (MF and RQ) belong 
to different lineages, thereby indicating the potential for 
reproductive isolation between these two adjacent populations 
(Qi et al, unpublished data). These genetic and distribution 
characteristics render P. forsythii an excellent candidacy for 
the assessment of ring species. Therefore, this study aimed to 
determine whether P. forsythii had a ring species-like divergence 
pattern by characterization of its genetic structure and gene 
flow. 


2. Materials and Methods 


2.1. Ethics statement All experiments were approved by and 
were in accordance with the Ethics Committee of the School of 
Life Sciences, Lanzhou University. 


22. Sampling and species distribution modelling This 
study included 210 specimens from 17 localities, encompassing 
majority of the distribution range of P. forsythii (Table S1 
and Figure 1). Field identification was based on external 
morphological characteristics (Zhu et al, 1999). 

The current potential of the species and suitable range were 
estimated using Maxent version 3.33 (Elith et al, 2011; Phillips et 
al, 2006). Field records and the literature yielded 29 P. fors ythii 
occurrences (Zhang et al, 2010) (Table S2). Environmental data 
included 19 bioclimatic variables from the WorldClim database 
(http://www.worldclim.org/; Hijmans et al., 2005) (Table S3); 
spatial resolution was 2.5 arc minute (-4 km). Model complexity 
was evaluated with "ENMeval" in R 3.61 (R Core Team 2017), 
through the determination of the Maxent-corrected AIC 
under different parameters. Thirteen variables were included 
in the final analysis. Species distribution modeling (SDM) was 
performed using Maxent, with occurrence data randomly 
segregated into 80% training and 20% testing data. This model 
used the default convergence threshold, 5,000 maximum 
iterations, and 10 replicates. An ROC curve close to 1 indicated 
a better model performance (Phillips et al., 2006). 
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2.3. Laboratory analyses DNA was extracted from muscle 
samples, using the TIANamp Genomic DNA Kit (Qiagen, 
Beijing, China) in accordance with the manufacturer's 
instructions. Three mitochondrial genes (ND2, ND4, and Cyt-b) 
were amplified. Primers were designed using PRIMER 6.0. PCR 
annealing temperatures are enlisted in Table S4. Sequences 
were assembled manually using SegMan 7.10 in DNA Star, 
aligned in BioEdit version 7.1.30 (Hall, 1999), and concatenated 
using Sequence Matrix 1.7.8 (Vaidya et al., 2011) for further 
analyses. All mitochondrial markers were deposited in the 
GenBank database (accession nos. MK803503- MK804072 and 
MN872429-MN872473). 

Genomic DNA was incubated at 37°C with Msel (New 
England Biolabs, NEB), T4 DNA ligase (NEB), ATP (NEB), and 
Meel Y adapter N containing barcode, and fragments of 350 
bp in size were isolated using a Gel Extraction Kit (Qiagen) in 
accordance with the manufacturer’s instructions. Thereafter, the 
genomic DNA of the samples was amplified via PCR, and the 
constructed libraries were sequenced using the Illumina HiSeq 
sequencer with a 150-bp paired-end sequencing strategy. The 
sequencing quality requirement for most bases was set at >Q20. 
Raw reads were filtered to eliminate reads containing adapters 
or reads of low quality, such that downstream analyses were 
then based on clean reads. The whole-genome sequencing data 
of Phrynocephalus vlangalii (GWHAAFC00000000) were used as 
the reference genome. Sequencing data were then aligned with 
the reference sequence using the BWA software (parameters: 
mem -t 4 -k 32 -M) (Li and Durbin, 2009). SAMTOOLS 
(mpileup -m 2 -F 0.002 -d 1000) was used to detect SNP of 
individuals (Li et al., 2009). 


2.4. Data analysis 

Phylogeographic analysis Bayesian analyses were performed 
using AIC in MrBayes 3.2.6 (Huelsenbeck and Ronquist, 
2001) with the GTR+I+G model selected by MrModeltest 
2.3 (Nylander, 2004). For the Bayesian analysis, 20 million 
generations with the Markov Chain Monte Carlo (MCMC) 
simulations were performed until the average SD of split 
frequencies decreased to <0.01; the initial 25% were discarded 
as burn-in frequencies. Simulations used the default model 
parameters starting from a random tree. Outgroups were 
P. axillaris and P. mystaceus. All tree files were visualized using 
Figtree 13.1. 

Population genetic analysis Pairwise Fs; (between- 
population differentiation) and isolation-by-distance were 
estimated using Arlequin 3.5 (Excoffier et al., 2005). In this 
program, genetic distance is determined as F,,/(1-F,,). Straight- 
line geographic distances and ring distances among the sites 
were used to analyze isolation by distance. Ring distances were 
determined as geographic distances, assuming that dispersal 
routes follow ring distribution rather than straight lines across 
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Figure 1 Geographical features and locations of sampling sites. For locality names, see Table S1. Dot color corresponds to popula- 
tion structure derived from the phylogenetic tree. Arrows represent direction of gene flow. All populations were divided into nine 
groups based on phylogeny: Groupl, MF; Group 2, HT, YT, YC, PS, CL; Group 3, YJS; Group 4, BC, XH, ALE; Group 5, LN, LT; 
Group 6, JS; Group 7, YL; Group 8, RQ; Group 9, HTR, CL315. 
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the uninhabitable desert. Phylogeographic analysis suggested a 
genetic "gap" between populations RQ and MF. Thus, MF was 
selected as the origin and RQ as the terminus to determine the 
ring distance among all populations (Qi et al, unpublished data). 
Correlations between genetic and geographic distances were 
determined using a Mantel test with 1000 permutations and 
Arlequin 3.5. 
Gene flow between populations Using maximum-likelihood 
migration rates, between-population gene flow was estimated 
using Migrate 3.7.2 (Beerli, 1997). This approach uses MCMC 
searches to account for unequally effective population sizes 
and asymmetrical gene flow (Beerli and Felsenstein, 1999). As in 
some lineages, populations were not continuously distributed 
geographically; rather, they were divided by populations 
from other lineages. The present grouping strategy was thus 
determined through the phylogenetic structure and geographic 
distribution of all populations. All populations were divided 
into nine groups, and gene flow was determined (Figure 1). 
Initial values of effective population sizes and gene flow were 
estimated from F;;. The analysis involved 10 short chains (500 
trees used out of 10,000 sampled) and three long chains (5000 
trees out of 100,000). Results from 10 runs were averaged and 
95% confidence intervals were summarized. 

Furthermore, population divergence and mixing patterns 
were inferred from P. forsythii SNP data in TreeMix 1.13 
(Pickrell and Pritchard, 2012; Wang et al., 2015). Custom Perl 
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scripts were used to enumerate the variation of all individuals 
into populations. Thereafter, a maximum-likelihood tree was 
estimated and migration events were superimposed over it, 
such that the tree with the added migration episodes accounts 
for a greater proportion of the overall genetic variance. We 
allowed for 10-11 migration events. The final results were 
visualized using R 3.6.1. 

To better estimate historic gene flow between populations 
from MF and RQ, we used IMa2 based on the IM model 
(Hey, 2010b) with mtDNA data. In this analysis, demographic 
parameters, including divergence time (t), effective population 
size of each extant population (q0 and q1) and the ancestor (q2), 
and migration rates (m0 and m1) between two groups, scaled 
by the mutation rate, were estimated. The maximum priors for 
the parameters used for IMa analyses were t = 5, q = 10, and m= 
100. In the simulations, we conducted each IMa simulation for 
100000 steps with a burn-in of 100000 steps. 


3. Results 


3.1. Prediction of species distribution SDM analysis revealed 
that this model was a good fit with an AUC of 209. The most 
suitable habitats for P. forsythii were around the Basin (Figure 
2), with the western region being more habitable than the 
eastern region. Finally, we observed a gap in the inhabitable 


region along the eastern region. 
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Figure 2 Distribution models for Phrynocephalus forsythii. Warmer colors represent greater likelihood of occurrence (red - 6, blue - 1). 
Scale bar indicates occurrence probability. The arrow points represent the distributed gap points. 
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3.2. Sequence data 

Mitochondrial DNA We obtained 210 mitochondrial gene 
sequences with an average length of 3538 bp and identified 
137 haplotypes (Table S1). The nucleotide composition was as 
follows: A = 0.351, C = 0.270, G = 0.132, and T = 0.247. 
Single-nucleotide polymorphisms Sequencing of the GBS 
library comprising 146 samples yielded 786522636 clean reads 
with an average of 5387141 reads per sample. We successfully 
mapped 755235460 reads, and the average mapping rate was 
96.05%. Finally, 544574 SNPs were filtered from 11710171 raw 
SNPs. 


3.3. Population genetic analysis 

Phylogeographic analyses Haplotypes derived from the 
Bayesian analysis of the 17 populations formed three distinct 
clusters: first, the population from Minfeng (the Minfeng 
Group), was separated from all other populations, and the 
other populations were grouped into two clades with non- 
overlapping distribution: six populations located to the north of 
the Tarim Basin formed a large distinctive group (the northern 
group: JS, LT, LN, YL, BC, XH, ALE, HTR, CL315 and RQ), and 
nine populations located in the south of Tarim Basin formed a 
second distinctive group (the southern group: HT, YC, YT, PS, 
CL and YJS) (Figure 1). 

Population analyses Pairwise F,r differed significantly 
from zero (P < 0.001) (Table S5). Pairwise F;; was the greatest 
between MF and RQ despite their adjacent distribution. 
Genetic distance and ring distances were positively correlated 
(Mantel’s r = 0.407, P < 0.001), consistent with a pattern of 
isolation by distance around the ring (Figure 3). Furthermore, 
the correlation between straight-line geographic distances and 
genetic distance was significant (Mantel’s r = 0.243, P =0.02) 
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(Figure 3). 

Gene flow Most adjacent groups exhibited gene flow except 
for groups 1 (MF) and 8 (RQ) (Figure 1; Table S6). Concurrent 
with the mtDNA results, predicted migration events revealed 
frequent gene exchange between adjacent populations (Figure 
1). Integrating Migrate and TreeMix data revealed a gene- 
flow model of ring species, thereby indicating that almost all 
adjacent populations displayed gene flow except for MF and 
RQ. Furthermore, historic gene flow analysis based on IMa2 
indicated no gene flow between MF and RQ. The migration 
rates (2NM) from MF to RQ and vice versa were 0.05062 and 
0.04556, respectively. However, gene flow was confirmed only 
when 2NM values were significant in an LLR test (P « 0.05 
or less) (Hey, 2010a; Nielsen and Wakeley, 2001). The present 
results showed that the migration rate (2NM) from MF to RQ 
and vice versa was not significant. Therefore, we considered 
that there was no gene flow between MF and RQ (Table 1). 


4. Discussion 


Our results show that P. forsythii conforms to the ring species 
model. First, the P. forsythii recorded distribution and SDM data 
indicated that it had a continuous circular distribution around 
the Tarim Basin. Second, populations of MF and RQ are linked 
by a chain of interbreeding populations wrapping around the 
Tarim Basin. Finally, these two terminal populations, MF and 
RQ, are genetically distinct. 

The distribution of P. forsythii populations along the ring 
structured in a manner similar to that of the ring species 
model (Alcaide et al, 2014). Although SDM implied the presence 
of a gap in the eastern region of the basin, the distribution 
between MF and RQ was continuous. Therefore, we consider 
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Figure 3 Genetic distance (calculated from mtDNA) increases with geographic distance around the ring. A: Correlation when "ring 
distances" are used. B: Correlation when straight geographic distances between sites are used. We assumed no direct gene flow be- 
tween MF and RQ or across the center of the ring. Ring distances approximate potential gene flow between populations around the 


Tarim Basin. 
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the P. forsythii distribution to be ring-shaped. Such a ring- 
shaped structure may result from the development of a 
dense hydrological network around the Tarim Basin and the 
expansion of the Taklimakan desert. Phylogeographic analysis 
revealed that after populations first differentiated from MF, 
the western populations (eg, JS and YJS) diverged from the 
northern and southern groups, respectively. Accordingly, we 
speculated that P. forsythii originated from the northern side of 
the Kunlun Mountains. Thereafter, the remaining populations 
dispersed northward and eastward along the edge of Tarim 
Basin. The divergence time of the northern group and the 
southern group coincide with the formation of rivers and oases 
around the Tarim Basin and the expansion of the Taklimakan 
(Qi et al, unpublished data). These geological features probably 
explain the ring-shaped distribution structure observed herein 
in P. forsythii. Furthermore, the Taklimakan Desert is a barrier 
similar to that observed with Phylloscopus trochiloides, a classic 
ring species in Asia, which is isolated by a large treeless gap 
(Irwin et al, 2001). We thus speculate that after MF initially 
differentiated, the remaining P. forsythii populations expanded 


Yue Ql et al. 
Phrynocephalus forsythii 
as a Potential Ring Species 


317 


along the basin. 

Our Migrate and TreeMix results provided the strongest 
evidence for a gene flow pattern that supported P. fors ythii as 
a ring species. In particular, we found no break in gene flow 
throughout the ring of populations, except between MF and 
RQ (Figures 1 and 4). This generally uninterrupted gene flow 
is due to short-distance dispersal in a continuously distributed 
species (Joseph et al, 2008). Our previous studies on P. forsythii 
population history indicate that they have not experienced 
rapid expansion (Qi et al, unpublished data). Therefore, we can 
speculate that, after initially differentiating from the northern 
and southern groups, the westernmost populations gradually 
spread eastward. Therefore, the short-distance dispersal was 
probably a gradual process among P. forsythii that eventually 
led to gene flow between adjacent populations (Slatkin et al, 
1993). 

Finally, our Fery, migration and IMa analyses among 
populations consistently displayed restricted gene flow among 
terminal populations MF and RQ, suggesting some degree of 
reproductive isolation (Figure 1 and Table S5). As expected in 


Table 1 Maximum-likelihood estimates (MLE) and 90% highest posterior density (HPD) intervals of demographic parameters from IMa 


multilocus analyses 


10 q0 42 2N0M0 2N1M1 
MLEs 4.915 1.64 1.359 1.987 0.05062 0.04556 
HPD95Lo 4.747 1.101 0.663 1.959 0 0 

HPD95Hi 4.997 1.999 1.999 1.999 0.1681 0.1411 


t0: Time since major lineage divergence; q0, 41, q2: Effective population size of MF, RQ, and ancestral population; 2N0M0: Population 
migration rate from the RQ to MF; 2NIM I: Population migration rate from MF to RQ. 
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Figure 4 Inferred maximum-likelihood trees of mixture events. Abbreviations are the same as in Table S1. Graphs inferred for all 
populations, allowing ten (A) and eleven (B) migration events. Migration arrows are colored in accordance with their weight. 
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a ring species, MF and RQ are associated with populations 
within the chain. Although these other populations exhibit 
gene flow, we could also detect gradual and continuous genetic 
variation among them. The finding of a markedly strong 
correlation between genetic distance and geographic distance 
(Figure 3) indicates that variation in genomic composition 
will increase as populations move further away from each 
other during species dispersal. Adjacent populations would 
therefore possess very few genetic differences, whereas distant 
populations will accumulate enough differences, resulting 
in reproductive isolation (Slatkin, 1993). Previous studies 
have similarly suggested that long-distance dispersal in a 
continuously distributed species results in isolation by distance 
(Irwin et al., 2005). 

Owing to the lack of a direct evidence of reproductive 
isolation, we cannot definitively claim that P. forsythii is a 
ring species (Irwin et al, 2001; Pereira et al, 2011). However, the 
ring-shaped diversification of P. forsythii is notable as a clear 
manifestation of speciation processes, and we can consider our 
study populations as intermediate transition species (Niemiller 
et al, 2008; Nosil, 2008). Environmental changes may interrupt 
the ring-species model by breaking these populations into 
two or more species that do not exchange genes (Gavrilets et 
al., 1998; Doebeli and Dieckmann, 2003). Global warming, a 
known threat to P. forsythii (Sinervo et al, 2018; Qi et al, 2019), 
may cause the obliteration of some populations and interrupt 
the chain. Thus, a major challenge for modern evolutionary 
biology is to develop strategies that can protect the balance 
between gene flow and isolation mechanisms among potential 
ring species such as P. forsythii. In addition to these long-term 
implications, our work also deepens our existing understanding 
of species formation. 
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Appendix 


Table S1 P. forsythii samples information in this study. 


mtDNA SNP 
Sampling site Sample size Sample code Haplotype code 
number haplotype number 


Cele (CL) 13,14,15,16,17,18,19, h17,h22,h26,h33,h42, 
20,21,22,23,24,25,26 h49,h56,h63,h72,h83,h84 


43,44,45,46,47,48,49, 
50,51,52,53,54,55,56 


Hetian (HT) h102,h103,h104,h105,h106 


76,77,78,79,80,81,82, h122,h123,h124,h125,h126,h127, 


Yecheng (YO) 4 4 4 - 83,84,85,86,87,88,89 h128,h129,h130,h131,h132,h133 


98,99,100,101,102, 
Bachu (BC) 14 14 12 14 103,104,105,106,107,108 h3,h4,h5,h6,h7 


117,118,119,120,121, 
122,124,125,126,127 


Luntai (LT) 12 11 8 12 h11,h12,h13,h14,h15 


, 139,140,141,142,143,144,145, 
Yuli (YL) 11 11 5 11 147,149,150,151,152,153,154 h21,h23,h24,h25,h27,h28 


175,176,177,178,179,180,181,182,183, h44,h45,h46,h47,h48, 
184,185,186,187,188,189,190,191,192 h50,h51,h52,h53,h54,h55 


205,206,207,208, h67,h68,h69,h70, 
209,210,211,212,213 h71,h73,h74,h75,h76, 


Cele 315 (CL 315) 14 14 


Table S2 Coordinate point information for SDM analysis of P. forsythii. 


Longitude (E) Latitude (N) 
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Table S3 19 bioclimatic variables. 


Bioclimatic variables Latitude (N) 


Z 
=> 
R 


Mean Diurnal Range (Mean of monthly (max temp - 


min temp)) dd 


wW 
4 ga 
o 
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w 
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Temperature Seasonality (standard deviation *100) 36.96 


w 
[A 
o 
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Min Temperature of Coldest Month 37.62 


c 
= 
o 
oo 


Mean Temperature of Wettest Quarter 37.38 


w 
m 
o 
junk 
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Mean Temperature of Warmest Quarter 39.25 


Bio. 12 Annual Precipitation 40.04 

Bio-13 Precipitation of Wettest Month, 465 
Bio_14 Precipitation of Driest Month 41.81 

Bio1S Precipitation Seasonality (Coefficient of Variation) 4132 
Bio_16 Precipitation of Wettest Quarter 41.08 

Bio-17 = Precipitation of Driest Quarter 5 BM 
Bio_18 Precipitation of Warmest Quarter 41.5 


Table S4 Pair primers used for amplification and sequencing of the mitogenome gene fragment of P. forsythii. 


Primer name Sequence(5' to 3) Annealing temperature 


5'- GTAAATAAACTGATGCCTGC -3' 


5'- AGACCACGGGATTGCTAT -3' 


5'- AGGCGTAGGGAAGAGAAT -3' 


Table $5 The F,, among populations of P. forsythii. 
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L 0.76740 0.73730 0.63323 0.77226 0.52157 0.00000 


LT 0.79353 0.74620 0.63792 0.79031 0.68850 0.61510 0.88408 0.00000 


BC 0.35737 0.17124 0.57877 0.16799 0.83825 0.79627 0.79758 0.82269 0.76071 0.00000 


YC 0.92951 0.92455 0.88124 0.94106 0.94635 0.93670 0.95069 0.93917 0.94694 0.94373 0.92211 0.00000 
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Theta Values of 2Ney for each recipient group 


n co 
3 
ge 
ga 
[0] 


(2Nep) Group 1 Group 2 Group 3 Group 4 Group 5 Group 6 Group 7 Group 8 Group 9 


Group 2 0.0017 396.48 185.81 215.72 12.95 608.82 9.25 2.26 1.08 
Group3 — 00098 — 3109 — 195076288 4645 7144 — 683 —— M9 175 
Group 4 0.00212 302.16 294.05 158.48 16.01 301.41 7.28 3.78 2.93 
Grup — 0054 — 432 — 2634 — 3917 — ILI] — MOB — 70 14 175 
Group 6 0.00242 468.81 426.43 134.6 72.55 25.11 5.48 0 2.97 
Gow7 — 0008 — 15 — 0 — 0 20 20 00 20 39 
Group 8 0.02006 0 10.69 0 0 2.15 2.35 11.42 ABI 
Note: Group 1: MF; Group 2: HT, YT, YC, PS, CL; Group 3: YJS; Group 4: BC, XH, ALE; Group 5: LN, LT; Group 6: JS; Group 7: YL; Group 8 


RQ; Group 9: HTR, CL315 


